ENFORCING THE NON-NEGATIVITY CONSTRAINT AND MAXIMUM 
PRINCIPLES FOR DIFFUSION WITH DECAY ON GENERAL 
COMPUTATIONAL GRIDS 
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Abstract. In this paper, we consider anisotropic diffusion with decay, which takes the form 
a(x)c(x) — div[D(x)grad[c(x)]] — /(x) with decay coefficient a(x) > 0, and diffusivity coefficient 
D(x) to be a second-order symmetric and positive definite tensor. It is well-known that this partic- 
ular equation is a second-order elliptic equation, and satisfies a maximum principle under certain 
regularity assumptions. However, the finite element implementation of the classical Galerkin for- 
mulation for both anisotropic and isotropic diffusion with decay does not respect the maximum 
principle. Put differently, the classical Galerkin formulation violates the discrete maximum princi- 
ple for diffusion with decay even on structured computational meshes. 

We first show that the numerical accuracy of the classical Galerkin formulation deteriorates 
dramatically with an increase in a for isotropic media and violates the discrete maximum prin- 
ciple. However, in the case of isotropic media, the extent of violation decreases with the mesh 
refinement. We then show that, in the case of anisotropic media, the classical Galerkin formulation 
for anisotropic diffusion with decay violates the discrete maximum principle even at lower values 
of decay coefficient and does not vanish with mesh refinement. We then present a methodology 
for enforcing maximum principles under the classical Galerkin formulation for anisotropic diffusion 
with decay on general computational grids using optimization techniques. Representative numer- 
ical results (which take into account anisotropy and heterogeneity) are presented to illustrate the 
performance of the proposed formulation. 



1. INTRODUCTION 

In this paper we consider heterogeneous anisotropic diffusion with decay, which takes the form: 
a(x)c(x) — div[D(x)grad[c(x)]] = /(x) with a{x.) > and D(x) is a symmetric and positive definite 
second-order tensor. This equation is a hnear second-order elhptic partial differential equation |21] . 
There are many important problems in Mathematical Physics which give rise to this equation |60j . 
Also, this equation arises in numerical methods and mathematical analysis of transient problems 
[35] . Some of these cases include: 
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(a) For certain gases, the diffusion process is accompanied by a decay of the molecules of the 
diffusing gas, and the decay is proportional to the concentration of the gas. Such a phenomenon 
can be modeled as a diffusion equation with decay. 

(b) For certain problems, the governing equation of diffusion in a moving domain can be trans- 
formed into a diffusion equation with decay. 

(c) Application of the method of horizontal lines to the transient diffusion equation (which is a 
linear parabolic partial differential equation) gives rise to a diffusion equation with decay. 

1.1. Maximum principles and discrete maximum principles. From the theory of partial 
differential equations, it is well-known that the diffusion equation with decay satisfies a maximum 
principle under appropriate regularity assumptions. In some cases (but not always) the physically 
important condition that the concentration is non-negative is a direct consequence of a maximum 
principle. It is important to note that the classical maximum principle for diffusion with decay is 
different from the classical maximum principle for pure diffusion equation (see Theorem \2.1\ and 
Remark \ 2. 5\ in this paper). 

It is imperative that predictive numerical simulations employ accurate and reliable discretization 
methods. The resulting discrete systems must inherit or mimic fundamental properties of contin- 
uous systems. The non-negative constraint and maximum principles are some of the essential 
properties of diffusion- type equations. However, it is well-known (and also discussed below) that 
many numerical formulations (including the popular ones) may not give non-negative solutions or 
satisfy maximum principles for these types of equations on general computational grids. Another 
point to note is that the satisfaction of maximum principles and the non-negative constraint by a 
numerical formulation will be altered by the presence of the decay term. (That is, the conditions 
under which a numerical formulation satisfies maximum principles and the non-negative constraint 
for pure diffusion can be different from those for diffusion with decay.) This leads us to discrete 
maximum principles. 

The discrete analogy of a maximum principle is commonly referred to as a discrete maximum 
principle (DMP). Some factors that affect discrete maximum principles are: numerical formulation, 
mesh size, element type, nature of the computational domain (e.g., presence/absence of holes), and 
properties of the medium ~ decay coefficient, diffusivity coefficient, anisotropy, and heterogeneity. 

1.2. Prior numerical works. Numerous numerical formulations have been developed for both 
isotropic and anisotropic diffusion equations. These formulations are based on finite difference 
methods |44t [27] , finite volume method |5H [TOt [15] , finite element method [HI ? mixed method 
[ini[56l[171ll8l[9l[7l[8lllll[l5|, discontinuous Galerkin method [S E] , spectral element method 
|3U] . and mimetic method |26 [ 1341 [T^ 136] . Most of these methods can be extended to diffusion with 
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decay. However, none of the aforementioned specific formulations satisfy maximum principles (for 
both pure diffusion equation, and diffusion with decay). 

Lately, there is a surge in research activity on enforcing maximum principles, especially for 
diffusion- type equations. However, these earlier works differ from the proposed formulation as they 
have one or more following limitations: 

(a) The studies did not consider anisotropy and heterogeneity. It should be noted that develop- 
ing numerical formulations that satisfy for isotropic diffusion is much easier than anisotropic 
diffusion, and there are practical solutions to satisfy maximum principles under the classical 
Galerkin formulation for homogeneous isotropic medium. These include: 

(i) Any one-dimensional mesh with linear elements satisfies maximum principles under the 
classical Galerkin formulation. 

(ii) Any mesh with acute-angled triangles or (even right-angled triangles) will satisfy max- 
imum principles. Under certain milder restriction, a Delaunay mesh will also satisfy 
maximum principles. Now, with advances in computational geometry, software packages 
are available which can produce Delaunay meshes for reasonably complex geometries. For 
example, CGAL [I], Qhull [11[2], Triangle [59]. 

(iii) A mesh with rectangular elements with some restrictions on the aspect ratio satisfies the 
discrete maximum principle [15]. In particular, a mesh with square elements satisfy the 
discrete maximum principle. 

None of the above conditions (except Condition (i)) ensure the satisfaction of maximum prin- 
ciples and the non-negative constraint if the diffusivity tensor is anisotropic. 

(b) The studies did not consider the effect of decay. The decay term affects the classical maximum 
principle of second-order elliptic partial differential equation (see Theorem 12. ip . Moreover, the 
decay terms alters the conditions under which a formulation satisfies maximum principles. 

(c) The studies did not consider general computational grids, but instead derived conditions on 
the mesh and on the properties of the medium. They limited their studies to structured grids 
(rectangular elements, acute-angled triangles). 

We now briefly discuss some of the important works on discrete maximum principles. The earlier 

works on discrete maximum principles are from the finite difference literature. Some of these notable 

works are [641 [T6] . It is important to note that these studies did not consider anisotropy, and general 

computational grids. In References |23 1 166 1 [67 1 165]. sufficient conditions are derived for higher-order 

elements to satisfy discrete maximum principles, but the studies are restricted to one-dimensional 

problems or isotropic diffusion. Ciarlet and Raviart [17] considered isotropic diffusion with decay 

under the classical Galerkin formulation. The main goal of Reference [T7] is to get restrictions on the 

mesh to satisfy maximum principles, and not a methodology that works on general computational 
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grids. Herrera and Valocchi |22j have employed flow-oriented derivatives to enforce the non-negative 
constraint. However, the methodology is limited in scope as it is restricted to rectangular grids, 
and a special form of the diffusivity tensor. References [52l [371 ESI ES] addressed (pure) anisotropic 
diffusion using finite volume techniques. All these papers are some variants of the idea proposed 
by LePotier [52], which is to choose the location of sampling points for the concentration in each 
cell in such a way to meet the non-negative constraint. The methodology (which is proposed 
under the finite volume method) cannot be easily modified to fit into the framework offered by 
the finite element method (at least, not in the present form presented in these references), and 
till date, there is no extension of this idea to the finite element method. Some notable works on 
discrete maximum principles and monotonicity are in the Multi-Point Flux Approximation (MPFA) 
literature |43 [ I50 [ [3T]. and these works considered logically rectangular grids, or derived restrictions 
on the mesh and medium properties. 

Liska and Shashkov [32] proposed a non-negative formulation for pure anisotropic diffusion equa- 
tion based on conservative finite difference methods |58j . Nakshatrala and Valocchi [36] have ex- 
tended the variational multiscale and lowest order Raviart-Thomas mixed formulations to produce 
non-negative solutions based on optimization techniques. Also, Reference [46 [ Appendix] discusses 
various conditions to satisfy the non- negative constraint. Another interesting work is by Burman 
and Ern [13] who have derived a nonlinear stabilized Galerkin formulation that satisfies a discrete 
maximum principle on general grids but they considered isotropic diffusion. Other recent works on 
discrete maximum principle include \32\ [62t [28| [29] , and all these works focused on getting restric- 
tions on computational meshes to satisfy maximum principles. As discussed in Reference [46] . the 
idea of getting restrictions on the mesh and medium properties in the case of anisotropic medium 
as the conditions are stringent, and in some cases a mesh may not even exist. This paper is an 
extension of the ideas presented in References [Ml [3S] . 

1.3. Main contributions of this paper. The main contribution of the paper is to present a ro- 
bust methodology for enforcing maximum principles and the non-negative constraint for anisotropic 
diffusion with decay. The methodology is applicable for general computational grids with low-order 
finite elements. We also derive a (theoretical) sufficient condition for uniform computational meshes 
under which the classical Galerkin formulation for diffusion with decay satisfies the maximum prin- 
ciple for one-dimensional problems. 

1.4. An outline and symbolic notation used in this paper. The remainder of this paper 
is organized as follows. In Section [21 we present governing equations for anisotropic diffusion 
with decay, and clearly outline the problem statement. In Section [S] we present a methodology 
for enforcing the non-negativity constraint and maximum principles for anisotropic diffusion with 
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decay on general computational grids. In Sections [H we illustrate the performance of the proposed 
formulation using representative numerical examples. Finally, conclusions are drawn in Section O 
The symbolic notation adopted in this paper is as follows. Throughout this paper repeated indices 
do not imply summation. (That is, we do not employ Einstien's summation convention.) We shall 
make a distinction between vectors in the continuum and finite element settings. Similarly, we make 
a distinction between second-order tensors in the continuum setting versus matrices in the context 
of the finite element method. The continuum vectors are denoted by lower case boldface normal 
letters, and the second-order tensors will be denoted using L^T[h]X blackboard font (for example, 
vector X and second-order tensor D). In the finite element context, we shall denote the vectors 
using lower case boldface italic letters, and the matrices are denoted using upper case boldface 
italic letters. For example, vector v and matrix K. Other notational conventions adopted in this 
paper are introduced as needed. 

2. GOVERNING EQUATIONS AND PROBLEM STATEMENT 

Let Q C M"""^ be a bounded open set, where "nd" denotes the number of spatial dimensions. 
The boundary is denoted by dQ, which is assumed to be piecewise smooth. A spatial point is 
denoted by x G fi. The gradient and divergence with respect to x are denoted by grad[-] and div[-], 
respectively. The concentration of a chemical species is denoted by c(x). The boundary is divided 
into two parts: r° and such that r° U T"^ = dQ and r° n F^ = 0. F° is that part of the 
boundary on which Dirichlet boundary condition is prescribed, and F^ is the part of the boundary 
on which Neumann boundary condition is prescribed. The unit outward normal to the boundary 
is denoted by n(x). The governing equations for anisotropic diffusion with decay can be written as 
follows: 

(la) a(x)c(x) - div[B(x)grad[c(x)]] = /(x) in Q 

(lb) c(x) = cP(x) onF° 

(Ic) n(x) • D(x)grad[c(x)] = tP(x) on 

where a(x) > is the decay coefficient, D(x) is the diffusivity tensor, /(x) is the volumetric 
source/sink, cP(x) is the prescribed concentration on the boundary, and tP(x) is the prescribed fiux 
on the boundary. The diffusivity tensor is symmetric, and assumed to be bounded and uniformly 
elliptic. That is, there exists two constants < .^i < < +oo such that 

(2) eiyV < y^B(x)y < e2yV Vx G 17 and Vy G M"'^ 

Equation ([1]) is a second-order elliptic partial differential equation, and from the theory of partial 

differential equations, it is known to satisfy the following maximum principle: 

5 



Theorem 2.1 (maximum principle). Let c(x) G C^(0) n C(0), and a(x) G C°(fi) wit/i a(x) > 0. 
/n addition, div[B(x)] exists and is bounded in fi. 7/a(x)c(x) — div[D(x)grad[c(x)]] > m i/ien 
c(x) satisfies the following equation: 

(3) min c(x) > min c~(x) 
where 

(4) c~(x) := min(c(x),0) 

A proof to the above theorem can be found in any standard books on partial differential equations 
(e.g., References [20 ^ \21 \ [53]). Few remarks about the above theorem and its implications are in 
order. 

Remark 2.2. If c{x.) satisfies a(x)c(x) — div[D(x)grad[c(x)]] < in il. (and the remaining condi- 
tions in Theorem \2. 1\ hold) then c(x) satisfies the following equation: 

(5) max c(x) < max c'^{^) 
where 

(6) c'^(x) := max(c(x), 0) 

Remark 2.3. If a(x) < then the equation pa|) is called Helmholtz equation. It should be noted 
that Helmholtz equation does not satisfy a maximum principle similar to Theorem \2.1\ and a coun- 
terexample can be found in Reference [53]. This implies that the condition Q!(x) > in Theorem 
\2.1\ cannot be relaxed. 



Remark 2.4. It should be noted that one can find in the literature maximum principles even when 
c(x) does not belong to C'^(yt) (and even when c(x) is only measurable, for example see Reference 
[61] ). A detailed discussion of such results is beyond the scope of this paper, and is not central to the 
development of the proposed numerical formulation. An interested reader on maximum principles 
under weaker conditions can refer to [Ml [Ml [ZD [ID] o,nd references therein. 

Remark 2.5. For the case of pure diffusion (i.e., a(x) = Oj we have the following maximum prin- 
ciple. Let c(x) G C^(r2)nC(fi), and div[D(x)] exists and bounded in Q. . //— div[D(x)grad[c(x)]] > 
in $7 then c(x) satisfies 

(7) min c(x) = min c(x) 

Remark 2.6. It is important to note the difference in the maximum principles for pure diffusion 
(which is given in Remark \2.5\) and diffusion with decay (which is given by Theorem \2.1\) . In 
the case of a{x) > (that is, diffusion with decay), the "non-negative minimum" occurs on the 



boundary, whereas in the case of a{x.) = (that is, pure diffusion) the maximum principle says 
that the minimum occurs on the boundary. 

2.1. Consequence of maximum principles. Maximum principles have important mathematical 
consequences in the study of partial differential equations and physical implications in modeling. 
Maximum principles are often employed in proving well-posedness (in particular, uniqueness of 
solution), and obtaining point- wise estimates. For example, for Poisson's equation (which is a 
second-order elliptic partial differential equation) the uniqueness of solution is a direct consequence 
of the maximum principle [l2]. To illustrate an important physical implication, let us apply the 
maximum principle outlined above to the transient diffusive system given by equation ([T]) . We shall 
assume that r° = dn (that is, we prescribe Dirichlet boundary conditions on the whole boundary). 
If /(x) > (i.e., we have volumetric source), and cP(x) > (i.e., we have non-negative prescribed 
Dirichlet boundary conditions on the whole boundary); then from the maximum principle it can 
be inferred that the quantity c(x) is non-negative in the whole domain. That is, 

(8) c(x) > yxen 

Now, the question is whether a given numerical formulation gives non-negative solutions if the 
prescribed data on the boundary is non-negative and the prescribed forcing function is a source. 
Also, whether a chosen numerical formulation gives solutions that are in accordance with maximum 
principles. This leads us to the problem statement and the approach taken in this paper. 

Remark 2.7. Under certain conditions (on the forcing function and boundary conditions), the 
non-negative constraint can be a special case of a maximum principle as shown above. However, it 
should be noted that, in general, the non-negative constraint can be an independent result, and need 
not be a consequence of any known maximum principle. For example, one can construct a simple 
problem in which the non-negative constraint is not a consequence of the maximum principle given 
in Theorem \2.1\ To wit, one can have a forcing function that is a source in some region and a sink 
in some other region of the domain. For this case, the conditions given in Theorem \2.1\ are not 
met, but still one may have the non-negative constraint on the concentration of the diffusant. 

2.2. Problem statement and our approach. The problem statement can be written as follows: 
develop a finite element formulation for anisotropic diffusion with decay that satisfies the non- 
negative constraint and maximum principles on general computational grids for low-order finite 
elements. 

The proposed methodology is based on the following key idea. We start with the finite ele- 
ment formulation of the classical Galerkin formulation, which has a variational statement. To this 
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variational statement, we augment the bounds on the nodal concentrations given by the maxi- 
mum principle. The resulting problem belongs to convex quadratic programming, and is solved 
by the active-set strategy. The proposed methodology works for all low-order finite elements (e.g., 
two-node linear element, three-node triangular element, four-node quadrilateral element, four-node 
tetrahedron element, and eight-node brick element) as nodal concentrations satisfying the maxi- 
mum principle ensure that the maximum principle is met throughout the computational domain. 
The proposed methodology, in general, does not work for high-order elements as illustrated in 
Figure [H 

3. WEAK FORMULATION AND DISCRETE MAXIMUM PRINCIPLE 

Herein, we employ the classical (single-field) Galerkin formulation. We shall define the following 
function spaces: 

(9a) V := {c(x) G H^{^1) \ c(x) = cP(x) on r°} 

(9b) Q := {u'(x) G i7^(Jl) I u'(x) = Oonr°} 

where H^{p,) is a standard Sobolev space [TO]. For weak solutions, we can relax the regularity 
requirement on the diffusivity tensor D(x). We shall assume that each component of D(x) is square 
integrable, which is equivalent to saying that 



(10) / tr[D(x)^]D)(x)] dO < +oo 

n 



where tr[-] is the standard trace operator [TJj used in continuum mechanics. The classical Galerkin 
formulation for anisotropic diffusion with decay ([T]) reads: Find c(x) G V such that 

(11) B{w]c) = L{w) Vu;(x) G Q 

where the bilinear form and linear functional are, respectively, defined as 

(12a) B{w;c):= I grad[u'(x)] • B(x)grad[c(x)] df] + / i(;(x)a(x)c(x) dfi 

Jo. Jn 

(12b) L{w) := [ w{x)f{x)dn+ [ u'(x)tP(x)dr 

It is well-known that the above weak form (jlip is equivalent to the following variational statement 

(13) minimize -I3{c; c) — L{c) 

c(x)eP 2 

It may not be possible, in general, to obtain analytical solutions for equations pip - ()12p especially 
for realistic problems with complex geometries. In such situations one may have to resort to 



numerical solutions. Herein we employ the Finite Element Method (FEM). Let the domain be 
decomposed into "iVeZe" non-overlapping open element subdomains. That is, 

Nele 

(14) n=[jn^ 

e=l 

where a superposed bar denotes the set closure. The boundary of Vt^ is denoted as 50*^ := Cl'^ — Vt^. 
For a non-negative integer m, P'"(0'^) denotes the linear vector space spanned by polynomials up 
to m-th order defined on the subdomain il*^. We shall define the following finite dimensional vector 
spaces of V and Q: 

(15a) := {c'^(x) G V \ c'^(x) G C7°(n), c'^(x)|^, G ¥^{n^),e = I,- ■ ■ ,iVe/e} 

(15b) Q!" := {u;'^(x) G Q \ u;'^(x) G C'^{n),w\jc)\^, G P'^(O^), e = 1, • • • , iVe/e} 

where k is a non-negative integer. A corresponding finite element formulation can be written as: 
Find c'^ix) G V'^ such that 

(16) B{w'';c'') = L{w^) Vu;'^(x) G Q'^ 

3.1. A methodology for enforcing the non-negative constraint and maximum principles. 

Before we present a methodology for enforcing the non-negative constraint and (discrete) maximum 
principles under the classical Galerkin formulation, we present some definitions and relevant results 
from numerical optimization. We shall use the symbols ^ and ^ to denote component- wise in- 
equalities for vectors. That is, for given any two (finite dimensional) vectors a and h 

(17) a ^ b means that Cj < hi \/i 

Similarly one can define the symbol ^. Let us denote the standard inner-product in Euclidean 
spaces by < • ; • > ■ A problem in quadratic programming takes the form 

(18a) minimize fo{x) := — < x; Qx > — < x; g > 

(18b) subject to Ax ^ b (inequality constraints) 

(18c) Cx = d (equality constraints) 

The above problem belongs to convex quadratic programming if Q is positive semidefinite (which 
makes the objective function fo{x) to be convex). As the name suggests, convex quadratic pro- 
gramming is a special case of convex optimization. For further details on convex optimization and 
associated numerical algorithms see References (SJ HOj BO] . 

We now return to the finite element implementation of the classical Galerkin formulation of 
anisotropic diffusion with decay. After finite element discretization, the discrete equations take the 
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form 



(19) Kc=f 

where K is a symmetric positive definite matrix, c is the vector containing nodal concentrations, 
and / is the load vector (arising from the forcing function). The corresponding minimization 
problem can be written as 

(20) minimize — < c; Kc > — < c; f > 

where ^^ndofs" denotes the number of degrees of freedom in the finite element mesh (which is equal 
to the total number of nodes minus the number of nodes at which a Dirichlet boundary condition is 
enforced). As shown in Figures [3l and [TOt the finite element solution based on equation ()19p produces 
unphysical negative concentrations even for simple problems. A formulation corresponding to (j20p 
that satisfies the maximum principle (given by Theorem 12.11 and Remark 1 2. 2 p can be written as 

(21a) minimize — < c; Kc > — < c; f > 

(21b) subject to Cminl ^ C ^ Cmaxl 

where 1 is a vector of size ndofs containing ones, and Cmin and Cmax are, respectively, given by 
(22a) 

Cmin • — min c (x) where c (x) — min{c(x),0} 
(22b) Cmax := maxc'^(x) where c^(x) = max{c(x),0} 

A corresponding formulation that satisfies the non-negative constraint can be obtained by setting 
Cmin = 0, and omitting the upper bound (which is equivalent to the condition Cmax = +oo); and 
can be written as follows: 

(23a) minimize - < c; Kc > — < c; f > 

indofs 2 



(23b) subject to ^ c 

where is a vector of size ndofs containing zeros. 

Remark 3.1. The constraint Cminl ^ c ^ Cmaxl can be rewritten in the standard form given by 
equation ()18bp as follows: 

C Cmaxl 
C ^ Cjninl 

The constraint ()23bp can be put in the standard form by rewriting it as: —c ^ 0. 
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Comparing with equation (jlSp . it is evident that the above problems (j2ip and (j23p belong to 
convex quadratic programming. The first-order optimality conditions (which are given by Karush- 
Kuhn- Tucker conditions) corresponding to equation ()23p take the following form: 

(25a) Kc = f + X 

(25b) ctO 
(25c) A ^ 

(25d) AjCj = {i = 1, ■ ■ ■ , ndofs) 

where A is a vector of Lagrange multipliers enforcing the constraint ()25bp . Similarly, one can write 
first-order optimality conditions for the optimization problem given by equation (j2ip . 



Remark 3.2. The above set of equations ()25p is not linear because of the inequality constraints 
(|25bp and (j25cp and complementary conditions ()25dp . 

4. REPRESENTATIVE NUMERICAL RESULTS 

In this section, we illustrate the performance of the proposed non-negative formulation for the 
anisotropic diffusion with decay using representative one- and two-dimensional problems. In all our 
numerical experiments we have employed the standard active-set strategy [40] to solve resulting 
convex quadratic programming problems. In all our numerical simulations we have taken the 
violated nodes under the classical Galerkin formulation as the initial active-set. This choice is 
motivated by the numerical studies reported by Nakshatrala and Valocchi [l6j in which it has been 
shown that the initial active-set based on the violated nodes from the underlying formulation, in 
most cases, takes fewer active-set strategy iterations (than, say, empty set as the initial guess). 

4.1. One-dimensional problem. Consider the following simple one-dimensional problem with 
homogeneous forcing function: 

d^c 

(26a) ac(x) - — = in J] := (0, 1) 

(26b) c(x = 0) = c(x = 1) = 1 

with a > 0. The analytical solution to the above problem is given by 

/ N 1 - exp(-Va) t n \ ^ exp(Va) - 1 i /- \ 

(27) c(x) = —— — exp(Vax) + —=- — exp(-Vax) 

exp(T/a) — exp(— ■y/aj exp(-y/aj — exp(— yaj 

In Figure [21 the analytical solution is plotted for various values of a. As one can see from the 
figure, sharp boundary layers exist for higher values of a. For obtaining numerical results, the com- 
putational domain is divided into four equal-sized elements. The numerical results obtained using 

the classical Galerkin formulation are shown in Figure O and the classical Galerkin formulation 
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clearly violates the discrete maximum principle for higher value of a. From this figure it is also 
observed that the larger the value of a the larger is the violation of the discrete maximum princi- 
ple. However, for one-dimensional problems, the violation of discrete maximum principle decreases 
with mesh refinement (which is not true, in general, in higher spatial dimensions especially when 
anisotropy dominates). 

We have solved the above one-dimensional problem using the proposed formulation, and have 
employed the active-set strategy to solve the resulting convex quadratic programming problem. 
We have taken a = 1000, and have employed the same computational mesh as discussed above. 
Figure H] illustrates how the active-set strategy performed at various iterations, and the active- 
set strategy converged in three iterations. In Figure [5l we have shown the performance of the 
^''clipping procedure" in which all the negative nodal concentrations from the Galerkin formulation 
are chopped off by setting them to zero. As discussed in the caption of the figure, the proposed 
methodology performs better than the clipping procedure. Also, the clipping procedure does not 
have a variational basis, and is considered as a "variational crime." 

In Figure El we plot the number of iterations taken by the active-set strategy with respect 
to number of (finite element) nodes. For one-dimensional problems, the violation of the discrete 
maximum principle decreases with mesh refinement, and eventually there will no violation of the 
discrete maximum principle. This can be seen in Figure [6] as the number of active-set strategy 
iterations is zero for (sufficiently) finer computational meshes for various values of a. In Figure [TJ 
we have shown the convergence of the proposed formulation with respect to mesh refinement for 
various values of a, and the proposed formulation performed well. 

We now derive sufficient conditions for uniform meshes for one-dimensional problems under the 
classical Galerkin formulation to satisfy the maximum principle. We shall use the following results 
from Matrix Analysis |55t l63j: Given Ax = b with b ^ 0, sufficient conditions to ensure that x ^ 
are 

(a) positive diagonal entries: An > 0, 

(b) non-positive off-diagonal entries: Aij < Vi 7^ j, and 

(c) strict diagonal dominance by rows: \Aii\ > J^j^^^i \ 

Remark 4.1. Note that the aforementioned sufficient conditions to ensure x ^ are quite strin- 
gent, and weaker (sufficient) conditions can be devised. For example, weaker sufficient conditions 
that ensures x y are: A is invertihle, and all the entries in A~^ are non-negative. 

Another sufficient condition that can he used requires that the matrix A to he an M-matrix, which 
is widely used in the numerical studies on flux and slope limiters |33j and iterative linear solvers 

|57| . An M-matrix is a non-singular matrix whose off-diagonal elements are non-positive and all 
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entries of the inverse are non-negative. Note that there are many equivalent definitions for an 
M-matrix [b7\ I63j . and the definition we just outlined is quite suitable for our discussion. 

Note that an M-matrix, by definition, has all the entries in its inverse to be non-negative. It 
can be shown that a matrix with positive diagonal entries, non-positive off- diagonal entries, and 
strict diagonal dominance by rows is an M-matrix ^5]. We have employed the sufficient conditions 
outlined just above this remark as they are easy to verify, and also suffice our purpose. 



We shall apply the above mathematical result to equation (|19p . which arises from the finite 
element discretization of diffusion with decay using the classical Galerkin formulation. The compu- 
tational domain is discretized using equal-sized two-node linear finite elements, and let h denotes 
the size of an element. Since the forcing function is assumed to be a source (that is, /(x) > 0), and 
the prescribed Dirichlet boundary conditions are non-negative; and we have / ^ 0. To get suffi- 
cient conditions for non-negative nodal concentration, we need to assess the entries of the "stiffness 
matrix" K. The entries of the stiffness matrix for an intermediate node (say i-th node) after the 
finite element discretization using two-node linear element under the classical Galerkin formulation 
take the following form: 



(28) 



ah 
~6" 



1 2 1 



Cj+l 



D r 



-1 2 



Ci-l 

Ci 
Ci+l 



where D denotes the diffusivity coefficient. Since a > 0, D > and h > 0; the conditions on 
positive diagonal entries and strict diagonal dominance are satisfied automatically. The condition 
on non-positive non-diagonal entries yields the following equation: 



(29) 



h < 




In Figure El we compare the above theoretical prediction with numerical simulations for various 
values of a, and the prediction is excellent. 



Remark 4.2. For the case of pure diffusion (that is, a = 0), equation (|29p implies that any 
uniform mesh using two-node linear finite elements satisfies the maximum principle. However, for 
diffusion with decay, there is a constraint on the mesh size h, which is proportional to D/a. That 
is, for a fixed D, the element size has to decrease with an increase in the decay coefficient to meet 
the maximum principle. This result highlights one of the main differences between diffusion with 
decay and pure diffusion under the classical Galerkin formulation. 
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4.2. Two-dimensional problem with isotropic medium. Consider the following two-dimensional 
problem with homogeneous forcing function: 

(30) Qc(x, y) - div[D grad[c(x, y)]] =0 in := (0, 1) x (0, 1) 

with D is assumed to be the identity tensor (i.e, isotropic medium) and a > 0. The geometry 
and boundary conditions for this two-dimensional problem are shown in Figure [9j The analytical 
solution is given by 

^ _ exp(Vax) + exp(Vay) 
exp(^) 

In this numerical study, we have taken a = 500. In Figure [TUl we show the numerical results 
obtained using the classical Galerkin formulation and the proposed formulation on a coarse compu- 
tational mesh. The classical Galerkin formulation violates the maximum principle, and the proposed 
formulation produces physically meaningful non-negative concentration even on the chosen coarse 
computational mesh. Since the diffusivity tensor is isotropic and the mesh is based on right-angled 
isosceles triangles, the violation of maximum principle under the Galerkin formulation (if it oc- 
curs) is due to the decay term. Moreover, the violation vanishes with sufficient mesh refinement. 
This fact is illustrated in Figure [TT] wherein we have employed a finer computational mesh, and 
the Galerkin formulation satisfies the maximum principle. However, it should be noted that the 
classical Galerkin formulation violates the maximum principle on fine unstructured computational 
meshes, which is illustrated in Figure [T2j 

In Figure [T3\ we plot the number of iterations taken by the active-set strategy with respect to 
number of (finite element) nodes. As discussed earlier, since the medium is isotropic, the violation 
of the discrete maximum principle again decreases with mesh refinement, and eventually there will 
no violation of the discrete maximum principle. This can be seen in Figure [13] as the number of 
active-set strategy iterations is zero for (sufficiently) finer computational meshes for various values 
of a. In Figure UM we perform numerical convergence studies of the proposed formulation for 
various values of a, and the proposed formulation performed well. 

4.3. Two-dimensional problems with anisotropic medium. Consider anisotropic diffusion 
in a bi-unit square plate = (0, 1) x (0, 1). The anisotropic diffusivity tensor is taken as follows: 

^_ / cos(^) sm{6) \ f h \ / cos{e) -sm{6) \ 
\ -sm{e) cos{9) J \ k2 J \ sm{9) cos{e) J 

with 9 = Tr/6, ki = 10000 and k2 = 1- The forcing function is taken to be zero (that is, /(x, y) = 0), 
and the decay coefficient is taken as a = 1. We prescribed Dirichlet boundary conditions on the 
whole boundary, and the prescribed concentrations are as follows: the left, right and top sides of the 
computational domain have a prescribed concentration of zero (that is, d'{x = 0, y) = cP(x = 1, y) = 
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(32) 



Table 1. Two-dimensional problem with anisotropic medium: Violation of maxi- 
mum principle with respect to mesh refinement using three- node triangular elements. 



mesh 


# of negative nodes 


% of nodes violated 


6x6 


6 


19.44 


12 X 12 


34 


27.78 


18 X 18 


90 


30.86 


21 X 21 


127 


31.74 


31 X 31 


301 


33.71 


41 X 41 


546 


34.56 


51 X 51 


854 


34.58 


101 X 101 


3500 


35.53 



cP(x, y = 1) = 0), and the bottom side of the computational domain has a prescribed concentration 
of cP(x, y = 0) = sin(7rx). The prescribed data in this problem meet all the conditions in Theorem 
12.11 and from the maximum principle we can infer the following: 

(33) < c(x) < 1 Vx G 

The problem is solved using two different computational meshes, and the numerical results are 
shown in Figure [15] (for three-node triangular mesh) and Figure [16] (for four-node quadrilateral 
mesh) . The amount of the violation of the maximum principle spatially for various mesh refinements 
is illustrated in Tables [1] and [2] In Figure [T71 we have shown the variation of minimum concentration 
with respect to mesh refinement. From these figures, it is evident that the negative concentration 
(which is the indication of the violation of the maximum principle) reached constant values for both 
three-node triangular and four-node quadrilateral meshes, and the violation existed irrespective of 
the mesh refinement. This is the main difference between the violation due to the decay term and 
the violation due to anisotropy. The violation of the maximum principle due to the decay term 
decreases with respect to mesh refinement, and eventually vanishes with mesh refinement. This 
fact is further illustrated in Figure [TE[ which shows the number of iterations taken by the active-set 
strategy for various values of a. 

4.4. Two-dimensional problem with a square hole. The computational domain is a bi-unit 
square plate Q := (0,1) x (0,1) with a square hole of dimension [4/9,5/9] x [4/9,5/9]. On the 
outer boundary we prescribe cP(x, y) = 0, and on the inner boundary we prescribe cP(x, y) = 2. 
The forcing function is taken to be zero (that is, /(x, y) = 0). The diffusivity tensor is same as the 
one employed in the previous subsection (see equation (|32p ). The computational mesh employed 
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Table 2. Two-dimensional problem with anisotropic medium: Violation of max- 
imum principle with respect to mesh refinement using four-node quadrilateral ele- 
ments. 



mesh 


# of negative nodes 


% of nodes violated 


6x6 


6 


25.00 


12 X 12 


36 


29.17 


18 X 18 


92 


31.17 


21 X 21 


127 


31.52 


31 X 31 


291 


32.98 


41 X 41 


530 


34.68 


51 X 51 


853 


35.22 


101 X 101 


3462 


35.44 



in this numerical simulation is shown in Figure [T9l Numerical results obtained using the Galerkin 
formulation and the proposed formulation are shown in Figure [20l and the proposed formulation 
performed well. 

4.5. Two-dimensional problem with heterogeneous anisotropic medium. This test prob- 
lem is similar to the one proposed in Reference [52], which addressed pure anisotropic diffusion 
equation. This test problem is considered as a good benchmark problem for testing numerical 
formulations for violation/satisfaction of discrete maximum principle. In this test problem, the 
diffusivity tensor is anisotropic and heterogeneous (that is, it varies spatially), and is given by 

y2 + ex2 -(l-e)xy\ 
-(l-e)xy x2 + ey2 J 

where e = 10~^. The domain is a bi-unit square plate: Q = (0, 1) x (0, 1). Homogeneous Dirichlet 
boundary conditions are prescribed on the entire boundary. The forcing function is taken to be 
/(x, y) = 1 if (x, y) G [3/8,5/8] x [3/8,5/8], and zero otherwise. Since the forcing function is non- 
negative, and homogeneous Dirichlet boundary conditions are prescribed on the whole boundary, 
from the maximum principle we have that the concentration is non-negative in the whole domain 
(that is, c(x) > in Cl). 

The numerical results that are obtained using the Galerkin formulation and the proposed for- 
mulation are shown in Figure EH The computational mesh that is employed in the numerical 
simulations is shown in Figure [22l Figure [23] shows the number of iterations taken by the active-set 

strategy under the proposed formulation for both three-node triangular and four-node quadrilateral 
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(34) 



Kx,y) 



meshes. In the case of heterogeneous anisotropic medium, violation of discrete maximum princi- 
ple occurs under the Galerkin formulation even for lower values of decay coefficient (e.g., even for 
a = 1), and the violation does not vanish even with mesh refinement. 

Remark 4.3. The proposed methodology works even for low-order three-dimensional finite elements 
like four-node tetrahedron element, eight-node brick element, and six-node wedge element. Herein, 
a three-dimensional problem is not solved as there are no computational challenges other than 
standard book-keeping. 

5. CONCLUSIONS 

In this paper, we have presented a methodology for enforcing the non-negative constraint and 
maximum principles for anisotropic diffusion with decay. The proposed method is obtained by 
adding constraints to the variational structure of the classical Galerkin formulation, and can handle 
general computational grids with low-order finite elements. The resulting equations form a convex 
quadratic programming problem, and are solved by employing the active-set strategy. Numerical 
experiments have shown that the rates of convergence with respect to mesh refinement in L^-norm 
and i?^-seminorm are about the same as for the original linear finite element method. Various 
representative numerical examples are presented to illustrate the good performance of the proposed 
formulation. 
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L2 element 



L3 element 




Q4 element qq g^g^ient 

Figure 1. This figure illustrates how the proposed methodology of enforcing the 
non-negative constraint and maximum principles works for low-order finite elements 
like two- node linear element (L2), three-node triangular element (T3), four-node 
quadrilateral element (Q4). The proposed methodology does not work for high-order 
elements like three-node quadratic element (L3), six- node triangular element (T6), 
nine- node quadrilateral element (Q9). In all the cases, the nodal concentrations are 
non-negative. For low-order elements, non-negative nodal concentrations ensures 
that the solution is non- negative within the whole finite element. In the case of 
high-order finite elements, enforcing non-negative nodal concentrations does not 
imply non-negative concentration throughout the element domain. 
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Figure 2. One-dimensional problem: Analytical solution for various values of alpha. 
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(c) a = 500, cmin = -0.1977 (d) a = 1000, Cmin = -0.2378 



Figure 3. One-dimensional problem: Comparison of the numerical solution ob- 
tained using the classical Galerkin formulation with the analytical solution for vari- 
ous values of the decay coefficient a. Note that the larger the value of a, the larger 
is the violation of the discrete maximum principle by the classical Galerkin formu- 
lation. 
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Figure 4. One-dimensional problem: This figure shows the variation of the numer- 
ical solution under the proposed formulation at various active-set strategy iterations 
for a = 1000. The active-set strategy converged after three iterations. Note that, 
for this problem, the converged numerical solution from the proposed formulation 
matches exactly at nodes with the analytical solution. 
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X 



Figure 5. One-dimensional problem: This figure compares the analytical solution 
with the numerical solution obtained using the "clipping procedure," which basically 
chops offs all the negative nodal concentrations obtained from the Galerkin formula- 
tion by setting them to zero. We have taken a = 1000 in this numerical simulation. 
The corresponding numerical solution obtained using the proposed methodolody 
is shown in Figure |4(d)| and the proposed methodolody performs better than the 
clipping procedure. 
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number of nodes number of nodes 

(c) a = 500 (d) a = 1000 



Figure 6. One-dimensional problem: These figures present the number of iterations 
required for the proposed formulation using the active-set strategy for various values 
of a with respect to the number of nodes. Note that the number of iterations 
required for the optimization to terminate increases as the value of a increases. 
After sufficient mesh refinement, there will be no violation of the discrete maximum 
principle, and there is no need to solve the constrained optimization problem. 
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Figure 7. One-dimensional problem: This figure presents numerical convergence 
of the proposed formulation with mesh refinement for various values of decay co- 
efficient. From the figure it is evident that the rates of convergence with respect 
to mesh refinement in L^-norm and i/^-seminorm are about the same as for the 
original linear finite element method. 
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Figure 8. One-dimensional problem: In this figure we compare the sufficient con- 
dition derived for uniform one-dimensional problems to satisfy maximum principles 
with numerical results, and the theoretical prediction is found out to be excellent. 
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Figure 9. Two-dimensional problem with isotropic medium: The forcing func- 
tion is taken to be zero. The analytical solution is given by c(x, y) = 
(1/ exp(-ya))(exp(-^/ax) -I- exp{^/a^y)). The Dirichlet boundary conditions are 
c(x, 0) = (1/ exp(^))(exp(^x) + 1) on Fi, c(l,y) = 1 + exp(^(y - 1)) on F2, 
c(x, 1) = exp(\/a(x — 1)) + 1 on F3, and c(0,y) = (1/ exp(-y/a))(l -|- exp(Y^y)) on 
F4. 
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(a) Three-node triangular mesh with 32 elements (b) Classical Galerkin formulation 
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(c) Proposed formulation 

Figure 10. Two-dimensional problem with isotropic medium: This figure shows 
the contours of concentration for a = 500 on a coarse computational mesh under 
the Galerkin formulation and the proposed formulation. Regions that have negative 
concentrations are indicated in white color. The proposed formulation produced 
physically meaningful non-negative concentrations in the entire computational do- 
main, Under the classical Galerkin formulation, approximately 24% of the total 
number of nodes have negative concentration. Also, under the classical Galerkin 
formulation, the minimum concentration is -0.4049, which occurred inside the do- 
main thereby violating the maximum principle. 
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(a) Three-node triangular mesh with 512 elements (b) Classical Galerkin formulation 

Figure 11. Two-dimensional problem with isotropic medium: This figure shows 
the contours of concentration for a = 500 on a fine computational mesh under the 
Galerkin formulation, and there is no violation of the maximum principle 
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(a) Unstructured mesh 



(b) Classical Galerkin formulation 
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(c) Analytical solution 



(d) Proposed formulation 



Figure 12. Two-dimensional problem with isotropic medium: The problem is 
solved on a fine unstructured mesh using the classical Galerkin formulation and 
the proposed formulation. Analytical solution is also shown in the Figure. The 
decay coefficient is taken to be a = 500. Regions that have negative concentrations 
are indicated in white color. The proposed formulation produced physically mean- 
ingful non-negative concentrations, and matched well with the analytical solution. 
Under the classical Galerkin formulation, approximately 14.2% of the total number 
of nodes have negative concentrations. Under the classical Galerkin formulation, the 
minimum value of concentration is —0.0466. Note that the negative concentration 
occurred mostly in the perturbed mesh region. 
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Figure 13. Two-dimensional problem with isotropic medium: These figures present 
the number of iterations required for the proposed formulation using active-set strat- 
egy at various values of a with respect to the number of nodes along each side of the 
computational domain (which is same in both x and y directions). Note that the 
number of iterations required for the active-set strategy to terminate increases as a 
increases. Again for this case, there is no violation of the discrete maximum prin- 
ciple after sufficient mesh refinement, and there is no need to solve the constrained 
optimization problem. 
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Figure 14. Two-dimensional problem with isotropic medium: This figure presents 
numerical convergence of the proposed formulation with mesh refinement for various 
values of decay coefficient. From the figure it is evident that the rates of convergence 
with respect to mesh refinement in L^-norm and i/^-seminorm are about the same 
as for the original linear finite element method. 
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Figure 15. Two-dimensional problem with anisotropic medium: The problem is 
solved using the Galerkin formulation (middle) and the proposed formulation (bot- 
tom). The left and right figures are, respectively, using 12 x 12 and 18 x 18 three-node 
triangular meshes. Regions that have negative concentrations are indicated in white 
color. Under the Galerkin formulation, 27.78% (for 12 x 12 mesh) and 30.86% of 
the total number of nodes have negative^ nodal concentration. The minimum con- 
centrations are —0.035 (for 12 x 12 mesh) and —0.022 (for 18 x 18 mesh). 
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Figure 16. Two-dimensional problem with anisotropic medium: The problem is 
solved using the Galerkin formulation (middle) and the proposed formulation (bot- 
tom). The left and right figures are, respectively, using 12 x 12 and 18 x 18 four-node 
quadrilateral meshes. Regions that have negative concentrations are indicated in 
white color. Under the Galerkin formulation, 29.17% (for 12 x 12 mesh) and 31.17% 
of the total number of nodes have nega^ve nodal concentration. The minimum 
concentrations are —0.020 (for 12 x 12 mesh) and —0.017 (for 18 x 18 mesh). 
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Figure 17. Two-dimensional problem with anisotropic medium: Variation of the 
minimum concentration with respect to mesh refinement for three-node triangular 
(T3) and four-node quadrilateral (Q4) meshes. 
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Figure 18. Two-dimensional problem with anisotropic medium: This figure shows 
the number of iterations taken by the active-set strategy under the proposed for- 
mulation for two different values of decay coefficient: a = 1 (top) and a = 500 
(bottom). The number of iterations are shown for both three-node triangular (left) 
and four-node quadrilateral (right) meshes. Equal number of nodes are employed 
along both x and y directions. Because of the anisotropy, the violation of the maxi- 
mum principle does not vanish with mesh refinement even for smaller values of decay 
coefficient (in this case, a = 1). 
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Figure 19. Two-dimensional problem with a square hole: Computational mesh 
using three- node triangular finite elements. 
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Figure 20. Two-dimensional problem with a square hole: Contours of the concen- 
tration obtained using the Galerkin formulation (left) and the proposed formulation 
(right) are shown in this figure. Regions that have negative concentrations are in- 
dicated in white color. The proposed formulation produced physically meaningful 
non-negative values for the concentration. Under the Galerkin formulation, approx- 
imately 26.92% of the total number of nodes have negative nodal concentrations. 
The minimum value of the concentration (which occurred inside the domain) is 
-0.0916. 
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Figure 21. Heterogeneous anisotropic medium: This figure shows the concentra- 
tion obtained using the Galerkin formulation (left) and the proposed formulation 
(right) for a decay coefficient of a = 1. The regions that have negative concentrations 
are indicated in white color. The proposed formulation produced physically mean- 
ingful non-negative values for the concentration. Under the Galerkin formulation, 
approximately 31.4% of the total number of nodes have negative nodal concentra- 
tions. The minimum value of the concentration is —0.0012. In the case of anisotropic 
medium, the violation of the maximum principle will occur even for smaller values 
of decay coefficient. Moreover, the violation, in general, will not vanish with the 
mesh refinement. 
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Figure 22. Heterogeneous anisotropic medium: Computational mesh using three- 
node triangular finite elements. 
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Figure 23. Heterogeneous anisotropic medium: This figure shows the number of 
iterations taken by the active-set strategy under the proposed formulation for two 
different values of decay coefficient: a = 1 (top) and a = 500 (bottom). The 
number of iterations are shown for both three-node triangular (left) and four-node 
quadrilateral (right) meshes. Equal number of nodes are employed along both x 
and y directions. Because of the anisotropy and heterogeneity, the violation of the 
maximum principle does not vanish with mesh refinement even for smaller values of 
decay coefficient (in this case, a = 1). 
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